home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / HyperCard Files / MacMathPascal / card_8077.txt < prev    next >
Encoding:
Text File  |  1988-09-05  |  4.4 KB  |  196 lines
  1. -- card: 8077 from stack: in
  2. -- bmap block id: 0
  3. -- flags: 0000
  4. -- background id: 2607
  5. -- name: SPLINE.PAS
  6.  
  7.  
  8. -- part contents for background part 22
  9. ----- text -----
  10. SPLINE.PAS
  11.  
  12. -- part contents for background part 23
  13. ----- text -----
  14. SPLINE.PAS
  15.  
  16. -- part contents for background part 17
  17. ----- text -----
  18. SPLINE(function)
  19.  
  20.  
  21.  
  22. PURPOSE 
  23.  
  24. The SPLINE routine smoothly interpolates a function, Y=Y(x).  The function is defined by a set of known values for Y at three or more values for the independent varable X. The  X values need not be uniformly spaced.  This function defines the natural cubic interpolatory spline.  In deriving the cubic interpolatory it is assumed Y is cubic in the interval between two points (Y'' is a constant).  It is further assumed that Y'' is zero at the end points of the function.  SPLINE works well for any class of functions for which the conditions are approximately true.  SPLINE uses the X, Y, and y'' values to perform the interpolation.  The y'' values can be found using the SPL procedure which can be found in the same file with SPLINE
  25.  
  26.  
  27. TYPE REQUIREMENTS
  28.  
  29. TYPE
  30.  FLOAT = EXTENDED or DOUBLE or REAL;
  31.  DATA = ARRAY[1..50] OF FLOAT;
  32.  FIXXED = INTEGER or LONGINT;
  33.  
  34. In addition, the procedure SPL is required.
  35.  
  36.  
  37. CALLING PROCEDURE
  38.  
  39. VAR
  40. X,F,S:DATA;
  41. Z:FLOAT;
  42. NDP:FIXXED;
  43. {X - INDEPENDENT VARIABLE ARRAY}
  44. {F - DEPENDENT VARIABLE ARRAY}
  45. {S - F'' ARRAY}
  46. {Z - INTERPOLATION POINT}
  47. {NDP - NUMBER OF DATA POINTS}
  48.  
  49.  
  50.  
  51.  
  52.  
  53.  
  54. SPL(X,F,S,NDP);{CALL SPL TO DETERMINE S}
  55. READLN(Z);{READ IN THE INTERPOLATION  POINT}
  56. Y:=SPLINE(Z,X,F,S,NDP);{FIND VALUE OF Y AT Z}
  57. WRITELN(Y);{PRINT OUT VALUE OF Y AT Z}
  58.  
  59.  
  60. EXAMPLE
  61.  
  62. The example program reads the number of data points (NDP).  It then reads in the data for two functions (X[1-NDP]), F1[1-NDP], and F2[1-NDP] and calls SPL for F1 and F2.  The program will then read in z values for which to find the interpolated F1 and F2 values and print the values until z=X[1].
  63.  
  64.  
  65. REFERENCES
  66.  
  67. Wendroff,B.: Theroretical Numerical Anaylsis, Academic Press, New York, 1966.
  68.  
  69.  
  70. -- part contents for background part 26
  71. ----- text -----
  72. 20
  73.  
  74. -- part contents for background part 6
  75. ----- text -----
  76. PROGRAM TEST_SPLINE;
  77. TYPE
  78. FLOAT = REAL;
  79. FIXXED = INTEGER;
  80. DATA = ARRAY[1..50] OF FLOAT;
  81. VAR
  82. I, IIN, NDP : FIXXED;
  83. ZIN : FLOAT;
  84. K : FLOAT;
  85. X2, F2, S2 : DATA;
  86. X : DATA;
  87. F : DATA;
  88. S : DATA;
  89. PROCEDURE see;
  90. VAR
  91. R : Rect;
  92. BEGIN
  93. HideAll;
  94. SetRect(R, 0, 35, 550, 330);
  95. SettextRect(R);
  96. Showtext;
  97. END;
  98. PROCEDURE SPL (VAR X, F, S : DATA;
  99. NDP : FIXXED);
  100. TYPE
  101. TEMP_STORE = ARRAY[1..65] OF FLOAT;
  102. VAR
  103. D : TEMP_STORE;
  104. E : TEMP_STORE;
  105. NM1, NM2, IP1, I, J, NSUB_INT : FIXXED;
  106. DOLD, DNEW : FLOAT;
  107. BEGIN
  108. NSUB_INT := NDP - 1;
  109. NM1 := NSUB_INT - 1;
  110. D[1] := X[2] - X[1];
  111. DOLD := (F[2] - F[1]) / D[1];
  112. FOR I := 1 TO NM1 DO
  113. BEGIN
  114. IP1 := I + 1;
  115. D[IP1] := X[I + 2] - X[IP1];
  116. E[I] := 2.0 * (D[IP1] + D[I]);
  117. DNEW := (F[I + 2] - F[IP1]) / D[IP1];
  118. S[IP1] := 6.0 * (DNEW - DOLD);
  119. DOLD := DNEW;
  120. END;
  121. S[2] := S[2] / E[1];
  122. NM2 := NM1 - 1;
  123. FOR I := 1 TO NM2 DO
  124. BEGIN
  125. IP1 := I + 1;
  126. D[IP1] := D[IP1] / E[I];
  127. E[IP1] := E[IP1] - SQR(D[IP1]) * E[I];
  128. S[I + 2] := (S[I + 2] - D[IP1] * E[I] * S[IP1]) / E[IP1];
  129. END;
  130. FOR I := 1 TO NM2 DO
  131. BEGIN
  132. J := NSUB_INT - I;
  133. S[J] := S[J] - D[J] * S[J + 1];
  134. END;
  135. S[1] := 0.0;
  136. S[NSUB_INT + 1] := 0.0;
  137. END;
  138. FUNCTION SPLINE (Z : FLOAT;
  139. X, F, S : DATA;
  140. NDP : FIXXED) : FLOAT;
  141. VAR
  142. ALPHA, BETA, HI, H1, S1 : FLOAT;
  143. I : FIXXED;
  144. FUNCTION QUBE (XIN : FLOAT) : FLOAT;
  145. BEGIN
  146. QUBE := XIN * XIN * XIN;
  147. END;
  148. BEGIN
  149. IF (X[1] < Z) AND (Z < X[NDP]) THEN
  150. BEGIN
  151. I := 1;
  152. WHILE Z > X[I] DO
  153. BEGIN
  154. I := I + 1;
  155. END;
  156. I := I - 1;
  157. HI := X[I + 1] - X[I];
  158. ALPHA := F[I + 1] / HI - S[I + 1] * HI / 6.0;
  159. BETA := F[I] / HI - S[I] * HI / 6.0;
  160. S1 := (S[I] * QUBE(X[I + 1] - Z) + S[I + 1] * QUBE(Z - X[I])) / (6.0 * HI);
  161. SPLINE := S1 + ALPHA * (Z - X[I]) + BETA * (X[I + 1] - Z);
  162. END;
  163. IF Z <= X[1] THEN
  164. BEGIN
  165. H1 := X[2] - X[1];
  166. S1 := (F[2] - F[1]) / H1 - H1 * S[2] / 6.0;
  167. SPLINE := F[1] + S1 * (Z - X[1]);
  168. END;
  169. IF Z >= X[NDP] THEN
  170. BEGIN
  171. H1 := X[NDP] - X[NDP - 1];
  172. S1 := (F[NDP] - F[NDP - 1]) / H1 - H1 * S[NDP - 1] / 6.0;
  173. SPLINE := F[NDP] + S1 * (Z - X[NDP]);
  174. END;
  175. END;
  176. BEGIN
  177. see;
  178. WRITELN(' HOW MANY DATA POINTS?');
  179. READLN(NDP);
  180. FOR IIN := 1 TO NDP DO
  181. BEGIN
  182. WRITELN(' ENTER X,F1,F2');
  183. READLN(X[IIN], F[IIN], F2[IIN]);
  184. END;
  185. SPL(X, F, S, NDP);
  186. SPL(X, F2, S2, NDP);
  187. ZIN := X[2];
  188. REPEAT
  189. WRITELN(' ENTER X');
  190. READLN(ZIN);
  191. K := SPLINE(ZIN, X, F, S, NDP);
  192. WRITELN('XIN =', ZIN : 15 : 7, ' Y =  ', K : 15 : 7);
  193. K := SPLINE(ZIN, X, F2, S2, NDP);
  194. WRITELN('XIN =', ZIN : 15 : 7, ' Y2 =', K : 15 : 7);
  195. UNTIL (ZIN = X[1]);
  196. END.